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Predicting  material  properties  of  disordered  systems  remains  a  long-standing  and  formidable  chal¬ 
lenge  in  rational  materials  design.  To  address  this  issue,  we  introduce  an  automated  software  frame¬ 
work  capable  of  modeling  partial  occupation  within  disordered  materials  using  a  high-throughput 
(HT)  first  principles  approach.  At  the  heart  of  the  approach  is  the  construction  of  supercells  con¬ 
taining  a  virtually  equivalent  stoichiometry  to  the  disordered  material.  All  unique  supercell  permu¬ 
tations  are  enumerated  and  material  properties  of  each  are  determined  via  HT  electronic  structure 
calculations.  In  accordance  with  a  canonical  ensemble  of  supercell  states,  the  framework  evaluates 
ensemble  average  properties  of  the  system  as  a  function  of  temperature.  As  proof  of  concept,  we 
examine  the  framework’s  final  calculated  properties  of  a  zinc  chalcogenide  (ZnSi-^Sea,),  a  wide-gap 
oxide  semiconductor  (Mg^Zni-^O),  and  an  iron  alloy  (Fei-^Cu^)  at  various  stoichiometries. 


INTRODUCTION 

Crystals  are  characterized  by  their  regular,  repeating 
structures.  Such  a  description  allows  us  to  reduce  our  fo¬ 
cus  from  the  macroscopic  material  to  a  microscopic  sub¬ 
set  of  unique  atoms  and  positions.  A  full  depiction  of 
material  properties,  including  mechanical,  electronic,  and 
magnetic  features,  follows  from  an  analysis  of  this  primi¬ 
tive  lattice.  First  principles  quantum  mechanical  calcula¬ 
tions  have  been  largely  successful  in  reproducing  ground 
state  properties  of  perfectly  ordered  crystals  [1,  2].  How¬ 
ever,  such  perfection  does  not  exist  in  nature.  Instead, 
crystals  display  a  degree  of  randomness,  or  disorder,  in 
their  lattices.  There  are  several  types  of  disorder;  includ¬ 
ing  topological,  spin,  substitutional,  and  vibrational  [3]. 
This  work  focuses  on  substitutional  disorder,  in  which 
crystallographically  equivalent  sites  of  a  crystal  are  not 
uniquely  or  fully  occupied.  Rather,  each  site  is  character¬ 
ized  by  a  statistical,  or  partial,  occupation.  Such  disorder 
is  intrinsic  in  many  technologically  significant  systems, 
including  those  used  in  fuel  cells  [4],  solar  cells  [5],  high- 
temperature  superconductors  [6,  7],  low  thermal  conduc¬ 
tivity  thermoelectrics  [8],  imaging  and  communications 
devices  [9],  as  well  as  promising  rare-earth  free  materials 
for  use  in  free  sensors,  actuators,  energy-harvesters,  and 
spintronic  devices  [10].  Hence,  a  comprehensive  compu¬ 
tational  study  of  substitutionally  disordered  materials  at 
the  atomic  scale  is  of  paramount  importance  for  optimiz¬ 
ing  key  physical  properties  of  materials  in  technological 
applications. 

Unfortunately,  structural  parameters  with  partial  oc¬ 
cupancy  cannot  be  used  directly  in  first  principles 
calculations — a  significant  hindrance  for  computational 
studies  of  disordered  systems.  Therefore,  additional  ef¬ 
forts  must  be  made  to  model  this  disorder.  One  method 
relies  on  the  reformulation  of  the  disordered  system  into 
an  average  effective  virtual  compound,  i.e.,  virtual  crys¬ 


tal  approximation  (VCA)  [11,  12].  In  the  VCA  approach, 
each  disordered  site  is  treated  as  a  virtual  atom  possess¬ 
ing  attributes  that  are  the  compositional-weighted  aver¬ 
age  of  the  actual  occupants.  An  advantage  of  the  VCA 
approach  is  that  the  computational  cost  of  a  disordered 
material  is  comparable  to  that  of  an  ordered  material. 
However,  this  approach  neglects  local  distortional  effects 
around  the  partially  occupied  sites — obscuring  fine  fea¬ 
tures  of  the  overall  structure. 

An  alternative  method  is  the  mean-field-type  coherent 
potential  approximation  (CPA)  [13],  often  implemented 
within  the  Korringa-Kohn-Rostoker  (KKR)  [14-16]  mul¬ 
tiple  scattering  formalism  for  enhanced  model  accu¬ 
racy  [17].  This  approach  simulates  interactions  between 
electrons  and  the  potential  via  propagators  (Green’s 
functions)  and  approximates  the  varying  potentials  of 
random  alloys  by  introducing  an  effective  medium  poten¬ 
tial  from  a  perfectly  ordered  lattice  [18].  Like  any  mean- 
field  approximation,  the  CPA  description  is  a  single-site 
approximation  and  thus  unable  to  resolve  short  range 
order  effects,  such  as  fine  resonance  structure  in  the  elec¬ 
tronic  density  of  states  (DOS). 

These  approaches,  among  others  [19-21],  grapple  with 
the  legitimacy  of  attributing  electronic  band  structure 
properties  to  random  alloys,  which  exhibit  no  transla¬ 
tional  long  range  order.  This  work  takes  a  different  ap¬ 
proach,  reformulating  the  issue  into  one  of  a  statistical 
nature.  While  not  all  material  properties  may  be  gar¬ 
nered  through  this  reformulation,  those  that  are  accessi¬ 
ble  are  of  significant  practical  importance,  including  the 
DOS,  band  gap  energy  E9Qp,  and  magnetic  moment  M. 

A  rigorous  statistical  treatment  of  substitutional  dis¬ 
order  at  the  atomic  scale  requires  utility  of  large  ordered 
supercells  containing  a  composition  consistent  with  the 
compound’s  stoichiometry  [22,  23].  However,  the  compu¬ 
tational  cost  of  such  large  supercell  calculations  has  tra¬ 
ditionally  inhibited  their  use.  Fortunately,  the  emergence 
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of  high-throughput  (HT)  computational  techniques  [24] 
coupled  with  the  exponential  growth  of  computational 
power  is  now  allowing  the  study  of  disordered  systems 
from  first  principles  [25]. 

Herein,  we  present  an  approach  to  perform  such  a 
treatment  working  within  the  HT  computational  frame¬ 
work  AFLOW  [26,  27].  We  highlight  three  novel  and  at¬ 
tractive  features  central  to  this  method:  complete  imple¬ 
mentation  into  an  automatic  high  throughput  framework 
(optimizing  speed  without  mitigating  accuracy),  utility 
of  a  novel  occupancy  optimization  algorithm,  and  use  of 
the  Universal  Force  Field  method  [28]  to  reduce  the  num¬ 
ber  of  DFT  calculations  needed  per  system. 


METHODOLOGY 

This  section  details  the  technicalities  of  representing  a 
partially  occupied  disordered  system  as  a  series  of  unique 
supercells.  Here  is  an  outline  of  the  approach: 

(1)  For  a  given  disordered  material,  optimize  its  partial 
occupancy  values  and  determine  the  size  of  the  derivative 
superlattices. 

(2)  (a)  Use  the  superlattice  size  n  to  generate  a  set  of 
unique  derivative  superlattices  and  corresponding  sets  of 
unique  supercells  with  the  required  stoichiometry.  (&) 
Import  these  non-equivalent  supercells  into  the  auto¬ 
matic  computational  framework  AFLOW  for  HT  first 
principles  electronic  structure  calculations. 

(3)  Obtain  and  use  the  relative  enthalpy  of  formation 
to  calculate  the  equilibrium  probability  of  each  supercell 
as  a  function  of  temperature  T  according  to  the  Boltz¬ 
mann  distribution. 

(4)  Determine  the  disordered  system’s  material  prop¬ 
erties  through  ensemble  averages  of  the  properties  cal¬ 
culated  for  each  supercell.  Specifically,  we  will  be  cal¬ 
culating  the  system’s  density  of  states  (DOS),  band  gap 
energy  Egap,  and  magnetic  moment  M. 

In  the  following  sections,  we  will  refer  to  a  model 
disordered  system,  Ag8  733Cd3  8Zr3  267,  to  illustrate  the 
technical  procedures  mentioned  above.  This  disordered 
system  has  two  partially  occupied  sites:  one  shared  be¬ 
tween  silver  and  zirconium,  and  another  shared  between 
cadmium  and  a  vacancy.  Working  within  the  AFLOW 
framework  [29],  we  have  designed  a  simple  structure  file 
for  partially  occupied  systems.  Adapted  from  VASP’s 
POSCAR  [30],  the  PARTCAR  contains  within  it  a  de¬ 
scription  of  lattice  parameters  and  a  list  of  site  coordi¬ 
nates  and  occupants,  along  with  a  concentration  toler¬ 
ance  (explained  in  the  next  section),  and  (partial)  occu¬ 
pancy  values  for  each  site.  To  see  more  details  about  this 
structure  or  its  PARTCAR,  please  see  the  Supplementary 
Materials. 


Table  I.  Evolution  of  the  algorithm  used  to  optimize  the  par¬ 
tial  occupancy  values  and  superlattice  size  for  the  disordered 
system  Ags.733Cd3.8Zr3. 267-  fi  indicates  the  iteration’s  choice 
fraction  for  each  partially  occupied  site,  (i  =  1,  2,  3,  . . . );  e; 
indicates  the  error  between  the  iteration’s  choice  fraction  and 
the  actual  partial  occupancy  value.  emax  is  the  maximum 
error  of  the  system. 
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1 

2 

1/2 

0.233 
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6 
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0 
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11 

12 

9/12 

0.017 
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12 

13 
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0.036 
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0.031 
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13 
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0.019 
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0.019 
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0.014 

0.019 

14 
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11/15 
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0 

0.00003 

15 

Determine  superlattice  size 

In  order  to  fully  account  for  the  partial  occupancy 
of  the  disordered  system,  we  would  need  to  generate  a 
set  of  superlattices  of  a  size  corresponding  to  the  low¬ 
est  common  denominator  of  the  fractional  partial  occu¬ 
pancy  values.  With  partial  occupancy  values  of  0.733 
(733/1000)  and  0.267  (267/1000)  in  the  disordered  sys¬ 
tem  Ag8.733Cd3.8Zr3.267,  we  would  need  to  construct  su¬ 
perlattices  of  size  1000.  Not  only  would  we  be  working 
with  correspondingly  large  supercells  (16,000  atoms  per 
supercell  in  our  example),  but  the  number  of  unique  su¬ 
percells  in  the  set  would  be  substantial.  This  would  ex¬ 
tend  well  beyond  the  capability  of  first  principles  calcula¬ 
tions,  and  thus,  is  not  practical.  It  is  therefore  necessary 
to  optimize  the  partial  occupancy  values  to  produce  an 
appropriate  superlattice  size. 

We  demonstrate  utility  of  an  efficient  algorithm  to  cal¬ 
culate  the  optimized  partial  occupancy  values  and  corre¬ 
sponding  superlattice  size  with  our  example  disordered 
system  Ag8.733Cd3  8Zr3.267  in  Table  I.  For  convenience, 
we  refer  to  the  algorithm’s  iteration  step  as  n '  ,  the  su¬ 
perlattice  index,  and  n  as  the  superlattice  size.  Quite 
simply,  the  algorithm  iterates,  increasing  the  superlat¬ 
tice  index  from  1  to  n '  until  the  optimized  partial  occu¬ 
pancy  values  reach  the  required  accuracy.  At  each  iter- 
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ation,  we  generate  a  fraction  for  each  partially  occupied 
site,  all  of  which  have  the  common  denominator  n '  . 
The  numerator  is  determined  to  be  the  integer  that  re¬ 
duces  the  overall  fraction’s  error  relative  to  the  actual 
site’s  fractional  partial  occupancy  value.  The  superlat¬ 
tice  size  corresponds  to  the  lowest  common  denominator 
of  the  irreducible  fractions  ( e.g .  see  iteration  step  8). 
The  maximum  error  among  all  of  the  sites  is  chosen  to 
be  the  accuracy  metric  for  the  system. 

For  the  disordered  system  Agg.733Cd3.gZr3.267,  given 
a  tolerance  of  0.01,  the  calculated  superlattice  size  is  15 
(240  atoms  per  supercell).  By  choosing  a  super  lattice 
with  a  nearly  equivalent  stoichiometry  as  the  disordered 
system,  we  have  reduced  our  supercell  size  by  over  a  fac¬ 
tor  of  60  and  entered  the  realm  of  feasibility  with  this 
calculation. 

Quite  expectedly,  we  see  that  the  errors  in  partial  oc¬ 
cupancy  values  calculated  for  silver  and  zirconium  are 
the  same,  as  they  share  the  same  site.  The  same  holds 
true  for  cadmium  and  its  vacant  counterpart  (not  shown) . 
Therefore,  the  algorithm  only  needs  to  determine  one 
choice  fraction  per  site,  instead  of  per  occupant  (as 
shown).  Such  an  approach  reduces  computational  costs 
by  guaranteeing  that  only  the  smallest  supercells  (both 
in  number  and  size)  with  the  lowest  tolerable  error  in 
composition  are  funneled  into  our  HT  first  principles  cal¬ 
culation  framework. 


Unique  supercells  generation 

With  the  optimal  superlattice  size  n,  the  unique  deriva¬ 
tive  superlattices  of  the  disordered  system  can  be  gener¬ 
ated  using  Hermite  Normal  Form  (HNF)  matrices  [31] . 
Each  HNF  matrix  generates  a  superlattice  of  a  size  cor¬ 
responding  to  its  determinant,  n.  There  exists  many 
HNF  matrices  with  the  same  determinant,  each  creat¬ 
ing  a  variant  superlattice.  For  each  unique  superlattice, 
we  generate  a  complete  set  of  possible  supercells  with 
the  required  stoichiometry  by  exploring  all  possible  oc¬ 
cupations  of  partially  occupied  sites.  However,  not  all 
of  these  combinations  are  unique — nominally  warrant¬ 
ing  an  involved  structure  comparison  analysis  that  be¬ 
comes  extremely  time  consuming  for  large  supercells  [31]. 
Instead,  we  identify  duplicates  by  estimating  the  total 
energy  of  each  supercell  in  a  HT  manner  based  on  the 
Universal  Force  Field  (UFF)  method  [28].  This  classical 
molecular  mechanics  force  field  approximates  the  energy 
of  a  structure  by  considering  its  composition,  connectiv¬ 
ity,  and  geometry,  for  which  parameters  have  been  tab¬ 
ulated.  Only  supercells  with  the  same  total  energy  are 
structurally  compared  and  potentially  treated  as  dupli¬ 
cate  structures  to  be  discarded,  if  necessary.  The  count 
of  duplicate  structures  determines  the  degeneracy  of  the 
structure.  Only  non-equivalent  supercells  are  imported 
into  the  automatic  computational  framework  AFLOW 


for  HT  quantum  mechanics. 

Supercell  equilibrium  probability  calculation 

The  unique  supercells  representing  a  partially  occu¬ 
pied  disordered  material  are  labeled  as  Si,  S2,  S3,  ..., 
S„.  Their  formation  enthalpies  (per  atom)  are  labeled 
as  Hp  \,  Hp 2,  Hp 3,  . . . ,  Hp4,  respectively.  The  forma¬ 
tion  enthalpy  of  each  supercell  is  automatically  calculated 
from  HT  first  principles  calculations  using  the  AFLOW 
framework  [26,  27].  We  take  the  supercell  with  the  lowest 
formation  enthalpy  as  a  reference  (ground  state  struc¬ 
ture),  and  denote  its  formation  enthalpy  as  Hp0.  The 
relative  formation  enthalpy  of  the  zth  supercell  is  cal¬ 
culated  as  A Hpi  =  Hp.i  —  Hpo  and  characterizes  its 
disorder  relative  to  the  ground  state.  The  probability 
Pi  of  the  zth  supercell  is  determined  by  the  Boltzmann 
factor: 

-AHFii/kBT 

□  _  i/*e 


i=l 

where  gt  is  the  degeneracy  of  the  zth  supercell,  AHF  l  is 
the  relative  formation  enthalpy  of  the  zth  supercell,  kp  is 
the  Boltzmann  constant,  and  T  is  a  virtual  “roughness” 
temperature.  T  is  not  a  true  temperature  per  se,  but  in¬ 
stead  a  parameter  describing  how  much  disorder  has  been 
statistically  explored  during  synthesis.  To  elaborate  fur¬ 
ther,  we  consider  two  extremes  in  the  ensemble  average 
(ignoring  structural  degeneracy):  1)  kpT  <  max  (A Hp^) 
neglecting  highly  disordered  structures  (AF^,;  0)  as 

T  — >  0,  and  2)  kpT  max  (A Hpi)  representing  the 
annealed  limit  (T  — >  00)  in  which  all  structures  are  con¬ 
sidered  equiprobable.  The  probability  Pi  describes  the 
weight  of  the  zth  supercell  among  the  thermodynami¬ 
cally  equivalent  states  of  the  disordered  material  at  equi¬ 
librium. 

Ensemble  average  density  of  states,  band  gap 
energy,  and  magnetic  moment 

With  the  calculated  material  properties  of  each  super¬ 
cell  and  its  equilibrium  probability  in  hand,  the  overall 
system  properties  can  be  determined  as  ensemble  aver¬ 
ages  of  the  properties  calculated  for  each  supercell.  This 
work  focuses  on  the  calculation  of  the  ensemble  average 
density  of  states  (DOS),  band  gap  energy  Egap,  and  mag¬ 
netic  moment  M.  The  DOS  of  the  zth  supercell  is  labeled 
as  1V,;(E)  and  indicates  the  number  of  electronic  states 
per  energy  interval.  The  ensemble  average  DOS  of  the 
system  is  then  determined  by  the  following  formula: 

n 

N (E)  =  J2Pi  x  W). 

2=1 
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Additionally,  a  band  gap  Egap ^  can  be  extracted  from 
the  DOS  of  each  supercell.  In  this  fashion,  an  ensemble 
average  band  gap  Egap  can  be  calculated  for  the  system. 
It  is  important  to  note  that  standard  density  functional 
theory  (DFT)  calculations  are  limited  to  a  description  of 
the  ground  state  [1,2].  Subsequently,  calculated  excited 
state  properties  may  contain  substantial  errors.  In  par¬ 
ticular,  DFT  tends  to  underestimate  the  band  gap  [32]. 
Despite  these  known  hindrances  in  the  theory,  we  demon¬ 
strate  in  the  next  section  that  our  framework  is  capable 
of  predicting  significant  trends  specific  to  the  disordered 
systems.  As  a  bonus,  the  calculation  of  these  results  are 
performed  in  a  high-throughput  fashion.  It  is  expected 
that  a  more  accurate,  fine-grained  description  of  the  elec¬ 
tronic  structure  in  such  systems  will  be  obtained  through 
a  combination  of  our  software  framework  and  more  ad¬ 
vanced  first  principles  approaches  [33-39]. 

In  the  same  spirit  as  the  N (E)  and  Egap,  we  consider 
the  calculation  of  the  ensemble  average  magnetic  moment 
M  of  the  system.  The  magnetic  moment  of  the  zth  super¬ 
cell  is  labeled  as  Mj.  If  the  ground  state  of  the  ith  struc¬ 
ture  is  non-spin-polarized,  then  its  magnetic  moment  is 
set  to  zero,  i.e.,  Mi  =  0.  Taking  into  account  the  impact 
of  signed  spins  on  the  ensemble  average,  this  approach 
is  limited  only  to  ferromagnetic  solutions.  Additionally, 
as  an  initialization  for  the  self-consistent  run,  we  assume 
the  same  ferromagnetic  alignment  among  all  of  the  spins 
in  the  system  (an  AFLOW  calculation  standard)  [26]. 
Finally,  the  ensemble  average  magnetic  moment  of  the 
system  is  calculated  with  the  following  formula: 

n 

M  =  ^Pix|Mi|. 

2  —  1 


EXAMPLE  APPLICATIONS 

To  illustrate  the  effectiveness  of  our  approach,  we  ana¬ 
lyze  disordered  systems  of  technological  importance:  zinc 
chalcogenides,  wide-gap  oxide  semiconductors,  and  iron 
alloys.  Unless  otherwise  stated,  the  supercells  used  in  our 
calculations  were  generated  with  the  lowest  superlattice 
size  nxct  needed  to  represent  the  composition  exactly. 


Zinc  chalcogenides 

Over  the  years,  zinc  chalcogenides  have  garnered  in¬ 
terest  for  a  dynamic  range  of  applications — beginning 
with  the  creation  of  the  first  blue-light  emitting  laser 
diodes  [43],  and  recently  have  been  studied  as  inorganic 
graphene  analogues  (IGAs)  with  potential  applications 
in  flexible  and  transparent  nanodevices  [44] .  These  wide- 
gap  II- VI  semiconductors  have  demonstrated  a  smoothly 


tunable  band  gap  energy  Egap  with  respect  to  compo¬ 
sition  [40-42],  Both  linear  and  quadratic  dependencies 
have  been  observed,  with  the  latter  phenomenon  referred 
to  as  optical  bowing  [45].  Specifically,  given  the  pseudo¬ 
ternary  system  A^Bi-^C, 

Egap(x)  =  [xeAc  +  (1  -  x)eBc }  -bx(l-x), 

with  b  characterizing  the  bowing.  While  Larach  et  al. 
reported  a  linear  dependence  ( b  =  0)  [40],  Ebina  et  al. 
[41]  and  El-Shazly  et  al.  [42]  reported  similar  bowing 
parameters  of  b  =  0.613  ±0.027  eV  and  b  =  0.457  ±0.044 
eV,  respectively,  averaged  over  the  two  observed  direct 
transitions. 

As  a  proof  of  concept,  we  utilized  our  developed  disor¬ 
dered  system  framework  to  calculate  the  compositional 
dependence  of  the  Egap  and  DOS  for  ZnSi-^Sea,  at  room 
temperature  (annealed  limit).  Overall,  this  system  shows 
relatively  low  disorder  (max(AFf^.j)  ~  0.005  eV),  ex¬ 
hibiting  negligible  variations  in  the  ensemble  average 
properties  at  higher  temperatures.  These  results,  illus¬ 
trated  in  Fig.  1,  are  directly  compared  against  experi¬ 
mental  measurements  [40-42] .  Common  among  all  three 
trends  (Fig.  1(a))  is  the  Egap  shrinkage  with  increasing 
xse >  ^  well  as  a  near  1  eV  tunable  Egap  range.  Our 
calculated  trend  demonstrates  a  non-zero  bowing  simi¬ 
lar  to  that  observed  by  both  Ebina  et  al.  [41]  and  El- 
Shazly  et  al.  [42].  A  fit  shows  a  bowing  parameter  of 
b  =  0.585  ±  0.078  eV,  lying  in  the  range  between  the  two 
experimental  bowing  parameters. 

We  also  plot  the  ensemble  average  DOS  plots  at  room 
temperature  for  x  =  0.00  (n  =  1),  0.33  (n  =  3),  0.67 
(n  =  3),  and  1.00  (n  =  1)  in  Figs.  l(b)-l(e).  The  plots 
echo  the  negatively  correlated  band  gap  relationship  illus¬ 
trated  in  Fig.  1(a),  highlighting  that  the  replacement  of 
sulfur  with  selenium  atoms  reduces  the  band  gap.  Specif¬ 
ically,  we  observe  two  phenomena  as  we  increase  the  con¬ 
centration  of  selenium:  (red  arrows)  the  reduction  of  the 
valence  band  width  (with  the  exception  of  xse  =  0.00 
(ZnS)  concentration),  and  (blue  arrows)  a  shift  of  the 
conduction  band  peak  back  towards  the  Fermi  energy. 
The  valence  band  of  ZnS  more  closely  resembles  that 
of  its  extreme  concentration  counterpart  at  xse  =  1-00 
(ZnSe)  than  the  others.  The  extreme  concentration  con¬ 
duction  peaks  appear  more  defined  than  their  intermedi¬ 
ate  concentration  counterparts,  which  is  likely  an  artifact 
of  the  ensemble  averaging  calculation. 

Finally,  we  consider  a  partial-DOS  analysis  in  both 
species  and  orbitals  (not  shown).  In  the  valence  band, 
sulfur  and  selenium  account  for  the  majority  of  the  states, 
in  agreement  with  their  relative  concentrations.  Mean¬ 
while,  zinc  accounts  for  the  majority  of  the  states  in  the 
conduction  band  at  all  concentrations.  Correspondingly, 
at  all  concentrations,  the  p-orbitals  make  up  the  majority 
of  the  valence  band,  whereas  the  conduction  band  con¬ 
sists  primarily  of  s-  and  p-orbitals.  These  observations 
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Figure  1.  Disordered  ZnSi-xSe^.  (a)  A  comparison  of  the  experimental  [40-42]  vs.  calculated  compositional  dependence  of 
the  band  gap  energy  Epa,p  at  room  temperature.  A  rigid  shift  in  the  E gap  axis  relative  to  experimental  results  of  ZnSe  (second 
ordinate  axis)  accounts  for  the  expected  systematic  deviation  in  DFT  calculations  [32].  Only  the  lowest  empirical  Effap  trends 
are  shown.  Error  bars  indicate  the  weighted  standard  deviation  of  the  ensemble  average  E gap.  (b)-(e)  Calculated  density  of 
states  plots  for  various  compositions:  (b)  x  =  0.00  (n  =  1),  (c)  x  =  0.33  (n  =  3),  (d)  x  =  0.67  ( n  =  3),  and  (e)  x  =  1.00 
(n  =  1).  The  straight  vertical  line  indicates  the  position  of  the  valence  band  maximum. 


are  consistent  with  conclusions  drawn  from  previous  opti¬ 
cal  reflectivity  measurements  that  optical  transitions  are 
possible  from  sulfur  or  selenium  valence  bands  to  zinc 
conduction  bands  [46]. 

Overall,  our  concentration-evolving  E90p  trend  and 
DOS  plots  support  a  continuing  line  of  work  [40-42] 
corroborating  that  this  system  is  of  the  amalgamation 
type  [47]  and  not  of  the  persistence  type  [46].  Notably, 
however,  reflectivity  spectra  shows  that  the  peak  posi¬ 
tion  in  the  E9Qp  for  ZnS  rich  alloys  may  remain  station¬ 
ary  [41],  which  may  have  manifested  itself  in  the  afore¬ 
mentioned  anomaly  observed  in  this  structure’s  valence 
band  width. 

Wide-gap  oxide  semiconductor  alloys 

Zinc  oxide  (ZnO)  has  proven  to  be  a  pervasive  material, 
with  far  reaching  applications  such  as  paints,  catalysts, 
pharmaceuticals  (sun  creams),  and  optoelectronics  [56]. 
It  has  long  been  investigated  for  its  electronic  proper¬ 
ties,  and  falls  into  the  class  of  transparent  conducting  ox¬ 
ides  [57].  Just  as  our  previous  zinc  chalcogenide  example, 
ZnO  is  a  wide-gap  II- VI  semiconductor  that  has  demon¬ 
strated  a  smoothly  tunable  band  gap  energy  E90p  with 
composition.  In  particular,  ZnO  has  been  engineered  to 
have  an  E90p  range  as  large  as  5  eV  by  synthesizing  it 
with  magnesium.  This  pairing  has  been  intensively  stud¬ 
ied  because  of  the  likeness  in  ionic  radius  between  zinc 
and  magnesium  which  results  in  mitigated  misfit  strain  in 
the  heterostructure  [58] .  While  the  solubility  of  MgO  and 


ZnO  is  small,  synthesis  has  been  made  possible  through¬ 
out  the  full  compositional  spectrum  [48-55] . 

As  another  proof  of  concept,  we  model  the  composi¬ 
tional  dependence  of  the  E9ap  and  DOS  for  Mga,Zni_xO 
at  room  temperature  (annealed  limit).  In  particular, 
we  chose  this  disordered  system  to  illustrate  the  breath 
of  materials  which  this  framework  can  model.  Similar 
to  ZnSi_xSea;,  this  system  shows  relatively  low  disorder 
(max (A Hp,i)  ~  0.007  eV),  exhibiting  negligible  varia¬ 
tions  in  the  ensemble  average  properties  at  higher  tem¬ 
peratures.  We  compare  our  results  to  that  observed  em¬ 
pirically  [48-55]  in  Fig.  2.  As  illustrated  in  Fig.  2(a), 
Ohtomo  et  al.  observed  a  composition  dependent  phase 
transition  from  a  wurtzite  to  a  rocksalt  structure  with 
increasing  XMg]  the  transition  occurring  around  the  mid 
concentrations.  We  mimic  this  transition  in  our  calcula¬ 
tions.  Empirically,  the  overall  trend  in  the  wurtzite  phase 
shows  a  negligible  bowing  in  the  Egap  trend,  contrasting 
the  significant  bowing  observed  in  the  rocksalt  phase.  We 
note  that  the  wurtzite  phase  E9Qp  trend  shows  a  slope  of 
2.160  ±  0.080  eV,  while  the  rocksalt  phase  shows  a  bow¬ 
ing  parameter  of  3.591  ±  0.856  eV.  Calculated  trends  are 
shown  in  Fig.  2(a).  Qualitatively,  we  also  observe  linear 
and  non-linear  Egap  trends  in  the  wurtzite  and  rocksalt 
phases,  respectively.  The  fits  were  as  follows:  we  observe 
a  slope  of  2.147  ±  0.030  eV  in  the  wurtzite  phase  and 
a  bowing  parameter  of  5.971  ±  1.835  eV  in  the  rocksalt 
phase.  These  trends  match  experiment  well  within  the 
margins  of  error.  We  observe  a  larger  margin  of  error 
in  the  rocksalt  phase,  particular  in  the  phase  separated 
region  (0.4  <  XMg  0.6).  This  may  be  indicative  of  the 
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Figure  2.  Disordered  Mg^Zni-^O.  (a)  A  comparison  of  the  experimental  [48-55]  vs.  calculated  compositional  dependence  of 
the  band  gap  energy  E3ap  at  room  temperature.  A  rigid  shift  in  the  Egrap  axis  relative  to  experimental  results  of  MgO  (second 
ordinate  axis)  accounts  for  the  expected  systematic  deviation  in  DFT  calculations  [32].  The  wurtzite  and  rocksalt  structures 
are  highlighted  in  blue  and  red,  respectively,  while  the  mixed  phase  structures  are  shown  in  black.  Error  bars  indicate  the 
weighted  standard  deviation  of  the  ensemble  average  E gap.  (b)-(e)  Calculated  density  of  states  plots  for  various  compositions: 
(b)  x  =  0.00  (n  =  1),  (c)  x  =  0.33  ( n  =  3),  (d)  x  =  0.67  ( n  =  3),  and  (e)  x  =  1.00  (n  =  1).  The  straight  vertical  line  indicates 
the  position  of  the  valence  band  maximum. 


significant  shear  strain  and  complex  nucleation  behavior 
characterizing  the  region  [49]. 

We  also  plot  the  ensemble  average  DOS  at  room  tem¬ 
perature  for  x  =  0.00  (n  =  1),  0.33  (n  =  3),  0.67  (n  =  3), 
and  1.00  (n  =  1)  in  Figs.  2(b)-2(e).  The  plots  not  only 
echo  the  positively  correlated  band  gap  relationship  il¬ 
lustrated  in  Fig.  2(a),  but  exhibit  the  aforementioned 
change  from  a  linear  to  non-linear  trend.  This  is  most 
easily  seen  by  observing  the  shift  in  the  conduction  band 
toward  the  Fermi  energy,  highlighted  by  the  blue  arrows. 
Contrasting  ZnSi_xSex,  we  do  not  observe  a  significant 
change  in  width  of  the  valence  band  as  we  vary  the  stoi¬ 
chiometry. 

Finally,  we  consider  a  partial-DOS  analysis  in  both 
species  and  orbitals  (not  shown).  Overall,  the  constant 
oxygen  backbone  plays  a  major  role  in  defining  the  shape 
of  both  the  valence  and  conduction  bands,  particularly  as 
XMg  increases.  This  resonates  with  the  strong  p-orbital 
presence  in  both  bands  throughout  all  concentrations. 
Zinc  and  its  d-orbitals  play  a  particularly  dominant  role 
in  the  valence  band  in  magnesium-poor  structures. 

Iron  alloys 

Despite  its  ubiquity,  iron  remains  at  the  focus  of  criti¬ 
cal  materials  research.  Even  as  new  phenomena  are  dis¬ 
covered  with  an  ever-growing  effort  to  explore  extreme 
conditions  [60-62],  there  exist  long-standing,  interest¬ 
ing  aspects  that  are  not  fully  resolved.  This  includes 


the  magnetic  character  of  the  (fee)  y-Fe  phase  at  low 
temperatures  [63-65],  among  other  complexities  in  its 
magnetic  phase  diagram  [66].  One  popular  approach  to 
studying  the  y-Fe  phase  is  through  the  Fei_xCux  dis¬ 
ordered  alloy  [63,  67,  68].  Nominally,  unary  copper 
and  iron  metals  with  fee  structures  are  nonmagnetic, 
but  together  exhibit  ferromagnetic  ordering  with  very 
high  magnetic  moments.  This  observation  has  led  to 
identification  of  Invar  and  anti-invar  behaviors,  which 
may  pave  the  way  to  enhanced  thermomechanical  actu¬ 
ators  [63,  67].  Fei_xCux  is  an  interesting  structure  in 
its  own  right,  as  it  has  extremely  low  miscibility  [69]. 
Overcoming  the  hurdle  of  developing  metastable  struc¬ 
tures  throughout  the  full  compositional  range  has  been 
the  focus  of  much  research  [70].  Such  metastable  struc¬ 
tures  have  demonstrated  novel  properties  like  high  ther¬ 
mal  and  electrical  conductivity  [71],  magnetoresistance, 
and  coercivity  [72]. 

As  a  final  proof  of  concept,  we  model  the  compositional 
dependence  of  the  magnetic  moment  M  for  Fei_xCux  at 
T  =  4.2  K  for  direct  comparison  against  experimental 
results  [59] .  Considering  both  the  sensitivity  of  magnetic 
properties  to  temperature  as  well  as  the  significant  dis¬ 
order  exhibited  in  this  system  (max (A Hp,i)  ~  1.63  eV), 
we  limit  our  analysis  to  the  low  temperature  limit.  This 
is  also  where  we  expect  our  framework  to  perform  op¬ 
timally,  which  considers  structures  relaxed  at  zero  tem¬ 
perature  and  pressure  [26].  The  results  are  illustrated 
in  Fig.  3.  Sumiyama  et  al.  shows  that  the  disordered 
system’s  phase  is  concentration  dependent,  with  a  phase 
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Figure  3.  Disordered  Fei-^Cua,.  (a)  A  comparison  of  the  experimental  [59]  vs.  calculated  compositional  dependence  of  the 
magnetic  moment  M.  Our  calculations  mimic  the  following  observed  phases  at  4.2  K:  x  <  0.42  BCC  phase,  0.42  <  x  <  0.58 
mixed  BCC-FCC  phases,  x  >  0.58  FCC  phase.  Error  bars  indicate  the  weighted  standard  deviation  of  the  ensemble  average 
E gap.  (b)  A  comparison  of  the  aforementioned  trends  with  calculations  performed  with  enhanced  superlattice  sizes  n.  (c) 
Calculated  unpolarized  density  of  states  (DOS)  plots  at  xcu  =  0.2  (n  =  5),  0.4  (n  =  5),  0.6  (n  =  5),  0.8  ( n  =  5). 


transition  from  bcc  to  fee  in  the  mid  concentrations  as 
xcu  increases.  Just  as  with  Mg^Zni-^O,  we  mimic  these 
observations  in  our  calculation.  The  overall  decreasing 
trend  in  M  with  reduced  Xpe  hr  Fig-  3(a)  matches  our 
expectations  well. 

With  such  a  simple  system,  we  also  explored  whether 
an  augmented  superlattice  size  n  enhances  our  results. 
While  the  concentration  remains  constant  for  n  above 
that  which  is  needed  for  the  desired  concentration  nxct, 
more  structures  are  introduced  into  the  ensemble  av¬ 
erage.  The  structures  themselves  also  increase  in  size 
by  a  factor  of  n  relative  to  their  parent  structure.  For 
xcu  =  0.2,  0.4,  we  simply  doubled  n  ( nenh  =  10),  while 
we  tripled  n  for  xcu  =  0.25  ( nenh  =  12).  With  only 
two  two-atom  structures  needed  to  describe  xcu  =  0.5 
at  n-xcti  we  were  able  to  increase  n  by  a  factor  of  five 
(nenh  =  10)  without  compromising  the  feasibility  of  the 
calculation.  A  comparison  of  results  calculated  at  nenh 
is  shown  in  Fig.  3(b).  At  most  concentrations,  we  ob¬ 
served  substantial  improvements  as  our  calculated  trend 
more  closely  follows  that  which  was  observed  empirically. 

Finally,  we  consider  this  system’s  ensemble  average 
DOS  in  Fig.  3(c).  In  general,  the  DOS  near  the  Fermi 
energy  decreases  with  increasing  Xcu,  with  some  insta¬ 
bility  near  the  mixed  phase  regions.  This  can  be  under¬ 
stood  using  the  Stoner  criterion  model  for  transitional 
metals  [73,  74].  Namely,  ferromagnetism  appears  when 
the  gain  in  exchange  energy  is  larger  than  the  loss  in 
kinetic  energy.  A  larger  DOS  at  the  Fermi  energy  in¬ 
duces  a  higher  exchange  energy  and  favors  a  split  into 
the  ferromagnetic  state.  The  competition  between  ferro¬ 


magnetic  and  paramagnetic  phases  can  be  inferred  from 
the  decreasing  M  trend  as  depicted  in  Fig.  3(a). 

ADDITIONAL  MATERIALS 

The  automated  approach  to  model  disordered  sys¬ 
tems  discussed  here  is  implemented  with  the  AFLOW 
framework  [26,  27]  and  can  be  downloaded  from 
http : //af lowlib . org. 

CONCLUSION 

In  this  work,  we  have  introduced  a  software  frame¬ 
work  capable  of  modeling  substitutionally  disordered  ma¬ 
terials.  Specifically,  the  framework  delivers  high  value 
properties  of  disordered  systems,  including  the  density 
of  states  (DOS),  band  gap  energy  Egap,  and  magnetic 
moment  M.  Though  a  number  of  technologically  signif¬ 
icant  examples,  we  have  illustrated  the  prowess  of  this 
highly  efficient  and  convenient  framework.  Such  materi¬ 
als  that  exhibit  highly  tunable  properties  are  of  critical 
importance  toward  the  goal  of  rational  materials  design. 
Without  loss  of  feasibility  or  accuracy,  the  framework  ex¬ 
ploits  highly  successful  high-throughput  first  principles 
approaches  in  more  complex,  real-world  systems. 
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